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Abstract 

A new diffuse interface model for a two-phase flow of two incompressible 
fluids with different densities is introduced using methods from rational con- 
tinuum mechanics. The model fulfills local and global dissipation inequalities 
and is frame indifferent. Moreover, it is generalized to situations with a sol- 
uble species. Using the method of matched asymptotic expansions we derive 
various sharp interface models in the limit when the interfacial thickness tends 
to zero. Depending on the scaling of the mobility in the diffusion equation we 
either derive classical sharp interface models or models where bulk or surface 
diffusion is possible in the limit. In the latter case a new term resulting from 
surface diffusion appears in the momentum balance at the interface. Finally, 
we show that all sharp interface models fulfill natural energy inequalities. 
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1 Introduction 

In recent years diffuse interface models have been successfully used to describe the 
flow of two or more immiscible fluids both for theoretical studies and numerical simu- 
lations. One fundamental advantage of these models is that they are able to describe 
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topological transitions like droplet coalescence or droplet break-up in a natural way. 
These ideas go even back to works of van der Waals and Korteweg, cf. [6] for a re- 
view and further references. In the case of two incompressible, viscous Newtonian 
fluids the basic diffuse interface model is the so-called "Model H" , cf. Hohenberg and 
Halperin [23]. It leads to the Navier-Stokes/Cahn-Hilliard system 

pd t v + p(v • V)v - div(2rj(c)Dv) + Vp = -be div(Vc ® Vc), (1.1) 

divv = 0, (1.2) 

d t c + v • Vc = div(mVp), (1.3) 

H = <7£~V'(c) - beAc. (1.4) 

Here p is the density, v is the mean velocity, Dv = |(Vv + Vv T ), p is the pressure, 
and c is an order parameter related to the concentration of the fluids (e.g. the 
concentration difference or the concentration of one component). Moreover, 77(c) > 
is the viscosity of the mixture, a is a constant related to the surface energy density, 
e > is a (small) parameter, which is related to the "thickness" of the interfacial 
region, ip is a homogeneous free energy density and p is the chemical potential. 
Capillary forces due to surface tension are modeled by an extra contribution eVc £g> 
Vc := eVc(Vc) T in the stress tensor leading to the term on the right-hand side of 
(11. II) . Moreover, we note that in the modeling diffusion of the fluid components is 
taken into account. Therefore div(mV/i) is appearing in (11.31) . where m = m(c) > 
is the mobility coefficient. 

One of the fundamental modelling assumptions is that the densities of both com- 
ponents as well as the density of the mixture p are constant. Of course, this restricts 
the applicability of the model to situations when density differences are negligible. 
Gurtin et al. [23] derived this model in the framework of rational continuum mechan- 
ics and showed that it satisfies the second law of thermodynamics in a mechanical 
version based on a local dissipation inequality. 

Lowengrub and Truskinovsky [2S] derived a thermodynamically consistent exten- 
sion of the Model H for the case of different densities, which leads to the system: 

pd t v + p(v • V)v - divS(c,£>v) + Vp = -de div(pVc <g> Vc), (1.5) 

d t p + div(pv) = 0, (1.6) 

pd t c + pv • Vc = div(m(c)Vp), (1.7) 

A* = -P' 2 tv + f V(c) - f div(pVc), (1.8) 

where S(c, Dv) = 2rj(c)Dv + A(c)divvl and A(c) is the bulk viscosity coefficient. 
Here the free energy has the density pa(e^ 1 ^/j(c)+e^^) per unit volume. A simplified 
version of this model has been successfully used for numerical studies, cf. Lee et 
al. [23 EH]- In contrast, there are - to the best of the authors' knowledge - no 
discrete schemes available which are based on the full model f ll.5p -( Tl~8|) . This may 
be due to fundamental new difficulties compared with Model H fll.ip - fll.4p . For 
instance, the velocity field v is no longer divergence free and therefore no solution 
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concept is available which avoids to determine the pressure p. At least analytically, 
these difficulties could be overcome, see Abels pQ for existence of weak solutions. 
Mathematically the coupling of the Navier-Stokes fll.5B - fll.6p and the Cahn-Hilliard 
part (ll.7p - fll.8p is much stronger, for instance the pressure p enters the equation for 
the chemical potential ( ll.Sp . and the linearized system is very different from the 
linearization of Model H, cf. Abels [2], where strong solutions locally in time are 
constructed. 

Alternative generalizations of the Model H for the case of different densities were 
presented and discussed by Boyer [8] and Ding et al. p2] . The model by Ding et al. 
consists of (jl.ip -( imp . but now for a variable density p = p(c). In order to justify 
this generalization they start from the mass balance equation 

d tPj + div(p jVj ) = (1.9) 

for the individual fluids j = 1,2 and define the mean velocity v of the mixture as 
volume averaged velocity v = «iVi + U2V2, where Uj is the volume fraction of fluid j. 
Then <^M> yields 

divv = 0, (1.10) 

cf. Section [2jbelow. In contrast to that Lowengrub and Truskinovsky define the mean 
velocity v as mass averaged/barycentric velocity pv = p x vi + P2V2, which yields 

d t p + div(pv) = 0. 

The incompressibility relation (11.101) of course has advantages with respect to nu- 
merical simulations - see the computations related to the model by Ding et al. in 

Em 

Unfortunately, neither global nor local energy inequalities are known to hold for 
(ll.ip - (ll.4p when p is not constant. The model by Boyer is more complicated. But 
it is derived using a volume averaged mean velocity, which leads to a divergence 
free mean velocity field too. The further derivation of Boyer differs from the one in 
[T2] and ours since the starting point are the equations for the conservation of linear 
momentum of each single fluid. Also for this model neither global nor local energy 
inequalities seem to be known, cf. also [7]. 

It is the purpose of the present paper to derive a thermodynamically consistent 
generalization of fll.lB - fll.4p to the case of non-matched densities based on a solenoidal 
velocity field v. More precisely, we will derive the system 

<9t(pv) + div(pv <g> v) + div(v <g> ^^m(<p)Vp) - div(2^(<p)Dv) + Vp 

= — as div(V(p ® V<p), 
div v = 0, 
d t (p + v ■ V<p = div(m(<p)Vp^), 

together with 
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where the order parameter (p = <p 2 — <fi stands for the difference of the volume 
fractions ipj, j = 1, 20. Here the free energy of the system has the density ae^f^) + 
ae^Y~ (per unit volume). We note that the first equation in the latter system is 
equivalent to 

pd t v + ((pv + Jj • Vj v - div(2n(ip)Dv) + Vp= -aediv(Vf <g> V<p), 

where J = pi ~ P2 m(y>)V/x, cf. Section [2] below. In comparison with the system derived 
in [12] there is the additional term (J • V)v in the equation for the linear momentum. 
This term vanishes in the case of matched densities, i.e., p = const.. But this term 
is crucial in the case of non-matched densities for consistency with thermodynamics. 
We note that in contrast to the model by Lowengrub and Truskinovsky (jl.5p - (jl.8p 
the usual continuity equation (11. 6p is not part of our system. Nevertheless there is 
conservation of mass in our system. More precisely, we have 

dtp + div (pv + j) = 0, (1.11) 

and in Section [2] it will be shown that in fact individual masses are conserved. Note 
that pj is the specific density of fluid j = 1,2 and that p = pl + p2 + P2 ~ Pl (p. We 
emphasize that according to equation (II. lip the (volume averaged) velocity v does 
not describe the flux of the density. In our model the flux of the density consists of 
the two parts: pv, describing the transport by the mean velocity, and a relative flux 
j _ — ^ 2 ~P 1 m (y})V/i related to diffusion of the components. Hence the diffusion of 
the components relative to the mean velocity leads to a diffusion of the mass density 
in the case that p~\ ^ p 2 . Moreover, we note that in the classical Model H effects 
related to diffusion of the components can play an important role and can lead to 
Ostwald ripening effects and disappearance of small droplets, cf. e.g. [34] . 

The structure of the paper is as follows: In Section [2] we will derive the generaliza- 
tion of Model H, described above in the framework of rational continuum mechanics. 
First we will use a local dissipation inequality and a choice of the energy flux as in 
[23] to derive restrictions for the form of the stress tensor and the chemical potential, 
which finally leads to our model after suitable constitutive assumptions. Then we 
briefly discuss the changes in the derivation if the energy flux is not specified at the 
beginning and Liu's Lagrange multiplier method is used, cf . [27] . In Section [31 we 
present a third approach to derive a thermodynamically consistent model, this time 
based on Onsager's variational principle. We consider the more general situation 
when either the system is subjected to gravitational forces or when one additional 
soluble species is present in both fluids. In the latter case, transport effects across 
the interface are taken into account, too, and we derive a diffuse interface analogue 
of Henry's law, cf. Section [31 

In Section H] we discuss the sharp interface asymptotics in the limit e — > for 
the diffuse interface model together with a soluble species. This is done by using the 

1 Other choices of order parameters are possible as well - see Section [5] 
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method of formally matched asymptotics. We show that the limit system depends 
essentially on the choice and the scaling of the mobility. Actually, we consider four 
cases related to choosing the mobility degenerate or non-degenerate and letting the 
mobility tend to zero or not. If the mobility m((p) vanishes as e — > 0, we end up with 
the classical model for a two-phase flow with the Young-Laplace law 

-[2r}Dv]±v + [p]j> = okv at T(t), 

where F(t) is the interface between the fluids, v is a unit normal to T(t), k is its mean 
curvature, a is a surface tension coefficient, and [.]! denotes the jump of a quantity 
at r(t) in the direction of v. Moreover, the interface is transported by the velocity 
of the fluid, i.e., 

V-i/-v = atr(t), 

where V is the normal velocity of T(t). If a soluble species with density w is present, 
we obtain the classical Henry condition for the jump of the concentrations of the 
soluble species. 

In the case of a constant mobility, we obtain in the limit e — > that the momentum 
balance in the bulk is given as 

fcdtv + div (v <g> (piV + ^^m V/i)) - ^Av + Vp = 0. 

In this case the diffusive flux — P2 2 Pl m V/i enters the momentum balance, cf. [5j. At 
the interface we obtain the Stefan type condition 

2(vi/-V) = mo\Vfi\t-u at 
for the evolution of the interface and again 

-[2r}Dv]±u + \p\±v = anu at T(t), 

for the jump of the stress tensor. Here \x satisfies 

2/i = ok — [iu]l at r(t), 

fi is harmonic in the bulk, and tuq > is a diffusion coefficient related to m. In 
particular the interface is no longer material and diffusion of mass through the bulk 
is still present in the model. In the case of a non-vanishing, degenerate mobility the 
evolution of the interface is governed by the surface diffusion law 

2(v • v - V) = mA r(i) /i, 

where Am) is the Laplace-Beltrami operator of T(t) and rh > is a diffusion coeffi- 
cient related to m. Moreover, the equation for the jump of the stress tensor is 

-^m((V r /i) • V r )v - [2r]Dv]tu + \p]±u = okv at T{t), 
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in this case. In particular the surface flux — P2 ~ pl mVrP enters the momentum balance 
at the interface. In Section [5] we prove that energy estimates are valid for sufficiently 
smooth solutions of the sharp interface models. Finally, several important identities 
for the formally matched asymptotics calculations are shown in the appendix. 

Acknowledgements: The authors are grateful to Hans Wilhelm Alt for several 
stimulating discussions and the telling us about [5], which resulted in the present 
model and helped to improve the presentation of the derivation of the model. More- 
over, we thank the anonymous referees for their valuable comments on the previous 
version of this paper. This work was supported by the SPP 1506 "Transport Pro- 
cesses at Fluidic Interfaces" of the German Science Foundation (DFG) through the 
grants GA 695/6-1 and GR 1693/5-1. 

2 Derivation of the Model 

In order to derive the diffuse interface model, one assumes a partial mixing of the 
macroscopically immiscible fluids in a thin interfacial region. Therefore let us intro- 
duce some terminology related to mixtures. In the following the fluids are labeled 
by j = 1,2 and they fill a domain Q C M. d . The total mass density of the mixture is 
denoted by p. Moreover, pj denotes the mass density of the fluid j, i.e., 

Mj = / Pj(x) dx 
Jv 

is the mass of the fluid j contained in a set V C f2 and we obtain p = p\ + p 2 - 
Moreover, we denote by Cj = ^ the mass concentration and note that C\ + c 2 = 1. 

Denoting by Jj the mass flux of fluid j, the mass balance equation in local form 
is given by 

d t pj + div Jj = 0. 

Defining the velocities Vj, j = 1, 2, of the single fluids as Vj = Jj/ pj the mass balance 
equation can be rewritten as 

d t pj + dw(pjVj) = 0. 

In what follows we assume that the volume occupied by a given amount of mass of 
the single fluids does not change after mixing, i.e., the excess volume due to mixing 
is zero. If pj is the specific (constant) density of the unmixed fluid j, we introduce 
Uj = j-. The assumption that the excess volume is zero results in 

«i + u 2 = l. (2.1) 
Expressed in terms of the mass concentrations C\ and c 2 , condition (12. ip reads as 
£iP + c 2 _p_ 1 i — £l + £2 

Pi P2 P Pi P2 
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Introducing the mass concentration difference c = c% — Cj, the above relation implies 
that p = p(c) with a function p defined via 



§(c+l) + ±(l-c) 



P(c) p 2 Pi 

We remark that possible choices for the order parameter in the phase field model 
are the mass concentration difference c, the density difference p := p 2 — p\ or the 
difference of volume fractions u := u<i — U\. 

We now introduce a suitable averaged velocity of the mixture. In contrast to the 
mass averaged/bary centric velocity v given by pv = piVi + p 2 v 2 , cf. Lowengrub and 
Truskinovsky [28], we choose the volume averaged velocity v of the mixture as in 
Boyer [8] and Ding et al. [12] . More precisely, we define 

Pi . P2 
V = M1V1 + W 2 V 2 = — Vi + — V 2 , 

Pi P2 

cf. [H [12]. As a consequence, we obtain, using the fact that the p/s are constants, 

div v = div(f vi) + div(f v a ) = div(| + |) = -ft (f + |) = -ftl = 0. (2.2) 

Furthermore, we denote by Jj = Jj — pjv the mass flux of the fluid j relative to the 
velocity v, i.e., 

d t pj + div(pj-v) + div = 0. 
Because of (12. 2p and pj = pcj = PjUj, we have 

dt(pcj) + v • V(pcj) + div Jj = 

or equivalently, 

d t Uj + v • Vuj + div 3j = 0, 

where Jj = Because of u\ +U2 — 1, we require 

Pj 

Ji + J 2 = ^ + ^ = 0, (2.3) 

Pi P2 

cf. (12.21) . In particular, we obtain 

d t (pc) + v ■ V(pc) + div J = 0, (2.4) 

where J = J 2 — Ji- In addition we have 

d t p = d t (pi + p 2 ) = - div(Ji + piv + J 2 + p 2 v) = - div(pv + J x + J 2 ). (2.5) 

If pi 7^ p 2 , we have in general div(Ji + J 2 ) 7^ 0. Hence the classical continuity 
equation does not hold with respect to the velocity v. This reflects the fact that we 
allow for mass diffusion in the system. 
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Instead of c, one could use u := u 2 — U\ as order-parameter. Because of 

p = pa + pc 2 = piu x + p 2 u 2 , (2.6) 
pc = pc 2 - pci = p 2 u 2 - pxut (2.7) 

and u 2 = u\ = we can also assume that p = p{u) and pc = pc{u). In order 
to have flexibility in the choice of the order parameter, we assume in the following 
that ip is any suitable order parameter such that p = p(f>) and pc = pc(<p) for some 
constitutive functions p, pc. Here pc = pc(<p) is the density difference. Possible 
choices are ip = Ui — u 2 as above, ip = c, and ip = pc. Then 

d t pc(ip) + v • Vpc(ip) + div J = (2.8) 

is equivalent to 

{d t cp + v • V<^) = - div J. (2.9) 

As in Gurtin et al. [23] , we assume that the inertia and kinetic energy due to the 
motion of the fluid relative to the gross motion is negligible. Therefore we consider 
the mixture as a single fluid with velocity v, which satisfies the law of conservation 
of linear momentum of continuum mechanics with respect to the volume averaged 
velocity. The density and the stress tensor are assumed to depend on additional 
internal variables like (p and Vy?. I.e., we assume that 

d t (pv) + div(pv ® v) = divf (2.10) 

for a tensor T, which has to be specified by constitutive assumptions. Here external 
forces are neglected for simplicity. 

The formulation (I2.10p can be rewritten as 

(dtp + v Vp)v + p(^v + (v • V)v) = divf. (2.11) 

In the above version the tensor T cannot be objective, i.e., an observer change 

(r,x*) = (t,a(t) + Q(t)x) 

will not lead to a transformation rule 

t = QT*Q T , 

see [S] for details. Using (12. 5p and div v = we obtain that the system f)2.1ip can be 
rewritten as 

p(d t y + (v • V)v) = div(T + v <g) (Ji + J 2 )) - ((Ji + J 2 ) ■ V)v . 

This system now allows for an objective tensor T = T + v <g> (J x + J 2 ). Using the 
fact that v is divergence free we introduce the unknown pressure p and restate the 
momentum equation as 

p(d t v + (v ■ V)v) = divS - Vp - ((Ji + J 2 ) ■ V)v (2.12) 

where S = T + pI = T + v® (J x + J 2 ) + pi, where I is the identity. 
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Remark 2.1 In the following we use the abbreviation 

J = Ji + J 2 

Because of (12. 3p . we have 

J = ^ij. (2.13) 
92 + Pi 

Then the total mass balance (12.51) can be rewritten as 

d t p + div(pv + J) = . (2.14) 

Using the mass balance we can restate the momentum balance (12.121) as 

pd t v+ ((pv + J) ■ V)v = divS-Vp (2.15) 

We note that pv + J = pv if v denotes the barycentric velocity defined by pv = 
Pi v i + P2 V 2- Moreover, ( I2.15P is equivalent to 

— / pv dx = [ TV ds x 



dt 



V(t) JdV{t) 



for any (sufficiently smooth) domain V(t) that is transported by v, i.e., the normal 
velocity of dV{t) coincides with u ■ v, cf (12.171) below. If V{t) is instead transported 
by the volume averaged velocity v, then (12. 15[) is equivalent to 



d_ 

dt 



/ pvdx= Tuds x - I (u-J)vds x = I TV ds x . (2.16) 

JV{t) JdV(t) JdV(t) JdVit) 



Finally, we assume the relative motion of the fluids to be diffusive, and we intro- 
duce a Helmholtz free energy density /(<p, V<p) (per unit volume). It will play the 
role of an interfacial energy for the diffuse interface. The total energy in a volume V 
is then obtained as the sum of the kinetic and the free energy, i.e., 

£y(<p,v)= / p( Lp )M-dx+ [ f(tp,V<p)dx= [ e(v,cp,Vcp)dx, 
Jv * Jv Jv 

where e(v, <p, V<p) := p(^)^- + /(<p, V<p). 

2.1 Derivation based on a Local Dissipation Inequality and 
Microstresses 

In the following V(t) C Q shall denote an arbitrary volume that is transported with 
the flow, i.e., the exterior normal velocity of V(t) is given by u ■ v(t) on dV (t). 
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In order to describe the change of the free energy due to diffusion, we introduce 
a chemical potential // and the outer unit normal v to dV(t) such that 

fi3-uds x — / ^-J-vds x 



2 

dV(t) JdV{t) 

is the energy transported into V(t) by diffusion and als x denotes integration with 
respect to the surface measure. Here the term — J dv ^ fiJ v ds x describes the change 

of free energy. The kinetic energy is transported by the velocity v + 3/p and the 
term — f gv , t \ 4j- J • v ds x describes the change of kinetic energy due to diffusion, see 
also [5], Section 10. 

Surface forces: Moreover, we assume the existence of a generalized (vectorial) surface 
force £ such that 

/ ip£-vds x 

JdV(t) 

represents the working due to microscopic stresses. Above and in the following ip is 
the material derivative d t ip + v ■ V<f. Finally, we note that 



(Tu) ■ v ds 



X 



dV(t) 



describes the working in a given volume V(t) due to the macroscopic stresses in the 
fluid. 

Second law of thermodynamics/local dissipation inequality: Similar as in |23j, we 
assume the following dissipation inequality, which is the appropriate formulation of 
the second law of thermodynamics in an isothermal situation: 



d_ 

di 



e(v, ip, Vip) dx 

V(t) 



< I (Ti/) ■ v ds x + I 'fk. v ds x — I /iJ • v ds x — J ^j- J ■ v ds x 

JdV(t) JdV(t) JdV(t) JdV(t) 

for every volume V(t) transported with the flow. This means that the change of 
total energy in time is bounded by the working due to macroscopic and microscopic 
stresses and the change of energy due to diffusion. 

We recall the transport theorem, see e.g. Liu [271 Theorem 2.1.]: 



2 ~ 



d_ 

dt 



[ fdx= f d t fdx+ f vvfds x = [ (d t f + div(vf))dx (2.17) 

JV(t) JV(t) JdV(t) JV(t) 



where V(t) is transported with the flow described by v. Therefore the equivalent 
local form is 

d t e + v • Ve - div(T T v) + div (^pj) - div(^) + div(//J) =: -V < 0. (2.18) 
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Using ffTTTUD . and <|27ID|> . we will simplify V. First of all, multiplying f[2TTUj) 

with v and using (I2.9p . we obtain 

a t (p^)+v-V(p^) = (^ P + vVp)^ + p(^v + (vV)v).v 

= -(divJ)^ + (divS- Vp-(J- V)v) -v 

= div (-±|v| 2 j + S T v-pv) -S : Vv. (2.19) 



Moreover, 



and 



4/ + v.V/=|£* + ^.(V„)- (2.20) 



<9pc 

div(<p£) — div(pj) = <p div £ + V(<p) • £ + Ptt~ <P — Vp • J 

<7<p 

= (|f p + div £)0 + (V^)' ■ £ + Vv : ( V<p <g> £) — Vp • J , 
because of (12.91) and where we have used 

d Xj (p = d t d Xj <p + v • Vd Xj <p + d Xj v ■ Vy? = (d Xj (p)' + d Xj v ■ V<p. (2.21) 
Thus we conclude that (12.181) is equivalent to 

+ (^^^ - ^ ■ (V<p)' - (S + V<p ® £) : Vv + Vp • J < 0, (2.22) 

where we have used div v = 0. 

In order to motivate the constitutive assumptions, we will derive some restric- 
tions for the constitutive relations specifying S,p, J,£ by an argument typical for 
rational continuum mechanics: To this end, we assume that S, J,£ are functions of 
Dy, tp, V<p, p, Vp only. Moreover, we assume that S is symmetric. Invoking general 
external forces and mass supplies in the equations, one argues that <p, <p, V<p, (V<p)", p, 
Vp, (Vv — | div vl) can attain arbitrary values for a given point in space and time 
and since / and £ do not depend on (p, (V<p)", we conclude from (I2.22p that 

^ Of .(<p,V<p) = 



necessarily. In particular, £ depends only on p, V<p. Hence (I2.22p reduces to 
_ fdf . dpc \ . 

-(S + V<p ® £) : Dv - (Vp <8>£) : ~(Vv- Vv T ) + V/i ■ J < 0, (2.23) 
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where Dv = |(Vv + Vv T ). Since the skew part of Vv can attain arbitrary values 
independent of Dv and 

(V? ® £) : -(Vv - Vv T ) = - (V95 <g> £ - £ ® V95) : -(Vv - Vv T ), 

we conclude that | (Vy? <S> $, — $, <E> Vy?) = and therefore 

|® Vy? = Vy9®£. 
Thus |£| 2 |Vy5| 2 = (Vv? • £) 2 and therefore 

df 

V^) = a(<p, V<p)V<p = g^{¥>, Vy?) (2.24) 

for some a((p, Vip). 

On the other hand, the first term after the equality sign in (I2.23P is linear in ip 
and /, J and /i are assumed to be independent of <p. Therefore the first term in (12.231) 
has to vanish to satisfy (I2.23P in general and 

dpc df ,. . . _ . 
—H = dw{a(cp, Vy?)V</?). 

Finally, the local dissipation inequality is satisfied if and only if 

V= (S + S/if®^- J : £>v- V/i- J > 0. (2.25) 
V dV<pJ 

Here S + Vy5 ® is also called viscous stress tensor since it corresponds to irre- 
versible changes of the energy due to friction in the fluids. 

Constitutive assumptions: Motivated by Newton's rheological law, cf. e.g. [23 Sec- 
tion 4.2.2], we assume that 

df 

S + Vif ® ^- = 27^)/^ 

for some function 77(95) > 0. 

Finally, we choose J (99, Vy5, /1, V/i) in the form 

3(<p,V(p,fi,Vfi) = -fh(<p)Vn, 

where m(tp) > 0, which corresponds to a generalized Fick's law. We remark that J 
can be chosen to be nonlinear with respect to V/i as long as f!2.25j) is fulfilled. 
Summing up, we derived the following diffuse interface model: 

pd t v + ((pv + J) ■ V)v - div (2i](<p)Dv) + Vp = - div(a(y>, Vy5)Vy5 <g> V<p), (2.26) 

divv = 0, (2.27) 
d t pc(ip) + v • V/5c(v?) = div(rh(if) V/i) (2.28) 
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together with 

^p = ^{ip, Vp) - div(a(p, V<p) Vp), (2.29) 

where a(<p, V<p) satisfies (12.241) and J = — ^~^ m(y)V/x due to f)2.13p . Assuming 
that the normal component of v and the normal derivative of p vanishes on the 
boundary of the fluid domain Q, J n pc(-,t) dx is a constant in time. By (12. 6p and 
(12.71) . we may express both mass densities as affine linear functions of pc. Hence, the 
total mass J n pi dx, i — 1, 2, of each liquid component is conserved. 

If <p = c is the mass concentration difference, then pc(c) = p(c)c, and we obtain 

pd t v + ((pv + J) • V)v - div (2r}(c)Dv) + Vp = - div(a(c, Vc)Vc ® Vc), (2.30) 

divv = 0, (2.31) 
d t {pc) + v • V(pc) = div(m(c)Vp) (2.32) 



together with 



9c 9c 
and J = -l^m(c)Vp. Here 

d(pb) 0-1 1 1 

— — = ap with a = — + — . 

oc 2pi 2p 2 

In the case that if = pc is the density difference, we have pc(<p) = <p, = 1, and 
therefore 

p<9 t v + ((pv + J) ■ V)v - div (2r)((p)Dv) + Vp = - div(a(<p, V<p) Vp <8> Vp), 

div v = 0, 
d t (p + v ■ V<p = div(m(<p)Vp) 

and 

df 

p = — (<p, V<p) - div(a(<p, V<p) Vp). 

Finally, in the case that tp = u 2 — U\ is the difference of volume fractions, we have 
pc(<p) = P2 ~ P1 + pi+pz (p anc i we obtain the system 

pd t v + ((pv + J) • V)v - div (2rj((p)Dv) + Vp = - div(a(<p, Vp)V<p <g> V<p), (2.34) 

divv = 0, (2.35) 
d t <p + v • Vp = — div J v , (2.36) 



together with 



df 

H v = -7±(<p, Vp) - div(a(<p, Vp)Vp) (2.37) 
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where we use a rescaled flux J v = { f 1 ^ 2 ) 1 J, and a rescaled chemical potential 
P<p = P1 + P2 ' fj,. We hence obtain 

J v = -m((p)V J = - £ ^= £1 m(v?)V/x ¥ , 

where m((p) = ( pi t p2 ) 2 rh((p). Usually we take the difference of volume fractions as 
order parameter. This has the advantage that the mass difference depends linearly on 
the order parameter and as usual in phase field models the values </? = ±1 correspond 
to unmixed "pure" phases. 

Finally we remark that (12.361) implies the mass balance 

d t p + div(pv + J) = 0. 

Using this we can rewrite the momentum balance (I2.34p as 

d t (pv) + div(pv <g> v) + div(v <g> J) - div(2r?(y?).Dv) + Vp = - div(a(y?, Vy?) V> <g> Vcp) . 

Remark 2.2 For the previous derivation it is essential to use the frame invariant 
version of the momentum equation (12. lip in order to derive a frame invariant model. 
If one does the derivation in the same manner using (I2.10p instead, the resulting 
model is 

dt(pv) + div(pv <g> v) - div (2rj{cp)Dv) + V> = - div(a(p, v» V> ® v», (2.38) 

divv = 0, (2.39) 
dtpc(tp) + v • Vpc(<p) = div(m(^)Vp) (2.40) 

together with 

cf. |3J for details. This model is not frame invariant since the chemical potential 

dp |v| 2 
dtp 2 



contains with the term jf-^r- a non-objective scalar 



2.2 Derivation based on a Local Dissipation Inequality and 
the Lagrange Multiplier Approach 

It is also possible to exploit a dissipation inequality without introducing generalized 
surface forces. We now assume a dissipation inequality 



— / e(v, (p, Vy?) dx + I J e • u ds x < 

v(t) ' ' Jav(t) 



dt 



for every volume V(t) transported with the flow for some general energy flux J e , 
which will be specified later. 
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Then the equivalent local form is 

d t e + v • Ve + div J e < 0. (2.42) 

Because of the conservation law (12.41) . we conclude that for every scalar function X v 
the inequality 

- V := d t e + v • Ve + div J e - X v {d t {pc) + v • V(pc) + div J) < (2.43) 
has to be valid. 

Using (12TTUD . (1239]) . f l2~2"Uj) we obtain that ( E^gj) is equivalent to 

+ ■ (Vy)' + div (-^J + T • v - A^J + J e ) - S : Vv + V\ v ■ J < 0, 

where we have used divv = 0. Making use of (12.211) we conclude that the latter 
inequality is equivalent to 

+ div^-^J + T-v-A v J + ^0 + J e ^j + VA^-J<0, (2.44) 

where again the equation div v = was used. 

If we assume now that S, A^,, and — ^-J + T-v — A (/3 J+ ^^0 — J e are independent 
of ip, we conclude that the first term after the equality sign in (I2.44p has to vanish 
for all values of (p. Thus 



d{pc) Of 

— =a dlv 

dip dip 



If we denote A^ = /i, then we obtain the same identity for the "chemical potential" 
fi as before. Moreover, if we now specify the energy flux as 

j e = M!j_T.v + "J-^, 
we end up with the local dissipation inequality 

V = (s + Vip <8> -^L^J : Dv — V/i • J > 0, 

which is just (I2.25p . Thus we can derive the model (I2.26p -( l2.29p as before by making 
the same constitutive assumptions. 
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3 Onsager's Variational Principle — a Third Ap- 
proach to Derive Thermo dynamically Consis- 
tent Models 

In this section, we follow a third pathway to obtain diffuse interface models in the 
spirit of fl2.34p - fl2.36l) and in agreement with the postulations of thermodynamics. It 
is based on Onsager's variational principle, see [30] and [13] . [19], [31] for applications 
in multi-phase flow. To widen the range of applications, we discuss two additional 
features. We allow for gravitational forces or, alternatively, we include the transport 
of a soluble species across fluidic interfaces as an additional effect. For simplicity, we 
refrain ourselves in the second case to a species that does not influence the surface 
tension at the interface. Adopting the notation from the previous sections, the order 
parameter will be denoted by (p and we take it to be the difference — U\ of the 
volume fractions Uj, j = 1, 2, of the two liquids involved. Hence, 

PVP) = — g — ^ + — 2 — ^ ' 

where p\ and p2 are the specific densities of liquid 1 and 2, respectively. By w, we 
denote the mass density of the soluble species which we assume to be dilute. 

Our starting point is the following set of general evolution equations, compare 
(EH and (J2II5D, 



P<9 t v + \^pv + ^ J • Vj v - div S + Vp = K, (3.2) 

div v = 0, (3.3) 
d t ip + v • Vip + div J v = 0, (3.4) 
d t w + v • Vw + div J w = (3.5) 

with p(-) as in IpTIJI and J v = J a rescaled flux, compare the discussion 

after f l2.37p . The stress tensor S is symmetric, and K denotes the force density. 
The additional equation (I3.5P is the mass balance of the soluble species and J w is 
the corresponding mass flux. These equations are supposed to hold in a space-time 
cylinder Q x (0, T) with Q C M. d being the domain where the process takes place. 
Conservation of mass requires that the normal components of and of 3 W vanish 
on dQ. As free energy, we choose 

f((p,V(p)dx+ / {g(w) + [3(ip)w}dx 
n Jn 

with f(z,p) = fi(z) + f2(p)- Here, g(-) is an entropic term, and /3((p) attains for 
ip < — 1 or </? > 1 the values f3\ or 02, respectively. In a sharp-interface limit, these 
parameters will reappear through the Henry jump condition 

?± = exp(/3 2 - fa) (3.6) 
w 2 
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at the interface separating the two phases (cf. Subsection 14.3 .3j) . The total energy is 
given by the sum of kinetic and free energy, hence 



,v| dx + / f(tp,Vtp)dx+ {g(w) + f3((p)w} dx. 
* Jn Jn 

The time derivative of the free energy is given as 

dF 1 f dp . l2 _ , f df _ , [ f df \ „ , 

-dt =2 y n a? 1 v| ^ " x + y n ^ " x - y n div j ^ rfx 

+ / p(<p)v ■ d t v dx + / (y'(iu) + /3(ip))d t wdx + / {3'((p)wd t ipdx, 
Jn Jn Jn 

where we assumed that the normal part of vanishes on dQ. Observe that 

d t v = ~ {k + div S - (Ypv + ■ Vj v - VpJ . 

Inserting the evolution equations and integrating by parts gives, assuming v = on 



d_ 

dt 



F = -\ I |^|v| 2 (v- Vy? + div J^dx- [ |^(v- Vy? + div J v )cfe 
^ Jn Of Jn o^P 

+ J div (J^ J (v • + div 3 V ) dx - J @'{ip)w (v ■ V<p + div J v ) dx 



dp 
dip' 



v • { K + div S - ( ( pv + J v ) • V ) v - Vp \ dx 



- I {g'{w) + /%)) (v ■ Vw + div 3 W ) dx 
Jn 

l v - v "{-| +d -(^)-^} & 

+ / v-Kcfe- / v • Vw {g'(w) + P(<p)) dx 
Jn Jn 



+ ■ Vpp dx — I Dv : S dx + / J w ■ V p w dx. 

Jn Jn Jn 

Here, we used 

V " ( ( PV + % J ^) " V ) V dX = ~\ i 1 V|2 (diV J " + VLf ' V) ^ 

which - as a consequence of the no-flux boundary conditions - follows easily by 
integration by parts. Moreover, we abbreviated 

to :=-div + 
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and 

fiu ■■= g'iw) + 0(<p) 

to denote the chemical potentials corresponding to ip and w, respectively. Recall 
that in general entropy production is due to external force fields and to gradients 
of velocity and of chemical potentials, see [TT], Chapter 3. However, if the specific 
densities of the external forces^] acting on the different species are identical, then 
those forces do not contribute to entropy production. Therefore, if no other external 
forces than gravity forces are applied, we may identify the rate of change of the 
mechanical work with 

dW f f f 

— — = — / v • Vipp^dx — / v • Vw// ffi dx + / v-K.dxdx. (3.7) 
dt Jn Jn Jn 

As 

J In Kgrav ■ v dx case with gravitational forces, 
dt 1 case with soluble species, 

equation ( 13. 7p is satisfied if 



K 



fipVcp + K grav case with gravitational forces, 
H v V(p + p w Vw case with soluble species. 



Here, K grav denotes the gravitational force. 

To determine the fluxes 3 V , 3 W and the stress tensor S, we introduce the dissipation 
functional 



2 |Q|2 |T |2 

VI 



^(3 V ,3 W ,S):= < , r / x +- L 7 L T+ ' \ > dx. (3.8) 
We use Onsager's variational principle which postulates 

This gives 3 V = -M(<p)Vfi v , 3 W = -K{<p)wVfi w , and S = 2q((p) Vv + VvT . Alto- 
gether, when a soluble species is around, we end up with 

dp 
dip' 

- div (2rj(ip)Dv) +Vp = /^V<p + n w Vw, (3.9) 

divv = 0, (3.10) 

d t ip + v • Vip - div (M(ip) V/i^) = 0, (3.11) 

d t w + v • Vw - div (K (^)wV^) = 0, (3.12) 

^=- div G^) + i + ^ (3 - 13) 

= g'(w) + p(<p) (3.14) 



P^v + ( ( pv + — J v 



2 i.e. the quotient of force density and mass density. 
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where Dv = |(Vv + Vv T ) and p{tp) = P2 2 pi (p + Pl + p2 . The evolution equation for 
w becomes 

9 ( w + v-Vw = div{K{<p)wV{g'{w) + /%))). (3.15) 
In the case of gravitational forces, the system has to be changed accordingly. 
Remarks: 

• If we choose the viscosity parameter in ( 13. 8 p as 7)((p, Dv), non-newtonian effects, 
e.g. shear thinning or shear thickening, can be included as well. 

• If we choose g(w) = w(\ogw — 1), the equations (l3.12p .( l3~7bT| result in the 
diffusion equation 

d t w + v ■ Vw = div(K((p)(Vw + u>V/%))). 

• The interfacial force term K = fi^Vip + fi w Vw can equivalently be written as 



V(f(<p, Vif) + g(w) + p(<p)w) - div Vip 



dVip 



Here, the first term is of pressure type whereas Vv?®^^ provides an additional 
stress tensor contribution representing interfacial forces. We can hence conclude 
that the derivations in Sections [2] and [3] up to a reinterpretation of the pressure 
lead to the same diffuse interface model if the dependence on w is omitted. The 
analogous observation for "Model H" has been already discussed in [23J. 

It is also possible to derive (13 . 9 p - ( |3TT4|) with the approaches discussed in Section 

m 

Taking the definitions of and K into account, we observe that K in fact 
does not depend on v. 



4 Sharp interface asymptotics 

In this section we identify the sharp interface limit of the diffuse interface model 
introduced in the preceding sections. We will use the method of formally matched 
asymptotic expansions where asymptotic expansions in bulk regions have to match 
with expansions in interfacial regions. There are four different asymptotic limits of 
interest. Two use a constant mobility and will either lead to a model where diffusion 
takes place through the bulk or to a model without any diffusion through the bulk, 
see also Abels and Roger [1] for the case when the densities in the two fluids are the 
same. Two cases are based on a mobility which is zero when the phase field takes 
the values ±1. In this case diffusion through the bulk is of course not observed in 
the sharp interface limit. But depending on the scaling we will either see surface 
diffusion along the interface or not. 
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4.1 The governing equations 

As usual for phase field models we introduce a scaling for / with respect to a small 
length scale parameter e as follows 

/(v?,Vv?) = ^|V^| 2 + -^) 
2 e 

where a is a constant related to the surface energy density. As in Section 3, we choose 
the difference of volume fractions as order parameter and we consider the following 
system 

dt(pv) + div(pv <g> v) + div(v <g> £^m(y?) V/i) - div(2r? (</?)£> v) + Vp (4.1) 

= — <7£div(V</? <8> V</?), 
divv = 0, (4.2) 
%> + (v • V)</? = div(m £ (ifi) V//) , (4.3) 

H=-ip'(<p)-aeA<p + p'(<p)w, (4.4) 
d t w + (v • V)iu = div(K(v?)«;V(^(w) + /%))). (4.5) 

To simplify the notation we drop the ip as index in the chemical potential. We assume 
that 

. p ( v ,) = ^a^ + a±£i, 

• g is convex, 

• /%) is smooth with 0(1) = /3 2 , = ft, 

• ?]((/?) is smooth and positive with 77(1) = r] 2 , = rji, 

• ip(ip) is a double- well potential such that = 1) = and ip(z) > if 
z 2 {1,-1}, 

• K(tp) is smooth, positive and such that K(l) = K 2 , K(—l) = K\. 
For the mobility m £ we distinguish four cases: 



m £ (ip) 



mo case I , 

em case II , 

^(1-V? 2 )+ caselll, 

,mi(l — y? 2 ) + case IV 



where mo, mi > are constants and (.)+ is the positive part of the quantity in the 
brackets. The total relevant energy in this scaling is 

F e (<p, v, w) = J |^|v| 2 + y|Vy?| 2 + -^W) + g(w) + 0(<p)w} dx . 
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For a solution (v e , p £ , ip e , p e , w £ ) of the system (I4.1l) - (l4.5p we perform formally matched 
asymptotic expansions. It will turn out that the phase field ip £ will change its values 
rapidly on a length scale proportional to e. For additional information on asymptotic 
expansions for phase field equations we refer to [HI [20] . 



4.2 Outer expansions 

We first expand the solution in outer regions away from the interface. We assume an 

expansion of the form v e = ^^ £ fc v fc , ip £ = YlkLti^Vk-, ^n expansion of (I4.4p 

in outer regions gives to leading order ip'{ipo) = and we obtain the stable solutions 
±1. We will denote by Q ± the regions where <po — ±1. In the cases II- IV the leading 
order expansion of the other equations are straightforward. We obtain: 

Pi (d t y + div(v g> v )) - 2 Vi div Dv + Vp Q = 0, (4.6) 

divv = 0, (4.7) 
d t w + (v • V)w = KiV ■ (w Vg'(w )) , (4.8) 

where i = 1,2 for x G and p\ = p(— 1) = pi,p2 = p(l) = p2- Due to the 

divergence free velocity we obtain 2 div Dv = Av and hence (14.61) simplifies to 

Pi(d t v + div(v <g> v )) - i]iAv + Vp = . 

We remark that (14. 8 p leads to the convection diffusion equation 

d t w + (v • V)w = KiAwo 

in the case that g(w ) = w(logw — 1). In the cases II-IV we will not need the chemical 
potential fi in the bulk. 

In case I the flux term — m(<p)V/U will enter the momentum balance (14. ip to 
leading order in the bulk. In addition also the evolution equation (14 .4p for the phase 
field will contribute to leading order to the sharp interface limit in the bulk. In fact 
in case I we obtain (14. 7p and (14 .8p together with 

Pi(d t v + div(v ® v)) +div(v® (^^)m Q V//) -^Av + Vp = 0, 

A/i = 0. 



4.3 Inner expansions 

We now make an expansion in an interfacial region where a transition between two 
phases takes place. 



4.3.1 New coordinates in the inner region 

We denote by T = (r(t)) t >o the smoothly evolving interface which we expect to be 
the limit of the zero level sets of ip when e tends to zero and will now introduce new 



21 



coordinates in a neighborhood of T. Choosing a time interval I cl and a spatial 
parameter domain U C IR d_1 we define a local parameterization 

7 : I x [/ -»■ R d 

of T. By we denote the unit normal to T(t) pointing into phase 2 (which is the 
phase related to ip — 1). Close to j(I x U) we consider the signed distance function 
d(t,x) of a point x to r°(t) with > if x E We now introduce a local 

parameterization of / x R d close to 7(1 x U) using the rescaled distance z = - as 
follows 

G £ (t, s, z) := (t, 7(t, s) + ezi/(t, s)) 
where s E U C R d_1 . We denote by 

V = 9 t 7 • 1/ 

the (scalar) normal velocity and observe that the inverse function (t, s, z)(t,x) := 
(G £ )~ l {t,x) fulfills 

To derive the last identity we used (2.6) and (2.20) of [10J. For a scalar function 
b(t, x) we obtain for b defined in the new coordinates via b(t, s(t, x), z(t, x)) = b(t, x) 
the identity 

^( ( ,x) = a,^6 + a lS .v,S + a,S = -ivai + h.o.t. 

where h.o.t. stands for higher order terms. With respect to the spatial variables we 
obtain, see Appendix, 

V x b = V v J+\d z bu (4.9) 

where Vr ez is the surface gradient on 

T £Z := {7(s) +ezu(s) \sEU} 

where here and in what follows we often omit the t-dependence. For a vector quantity 
x) written in the new coordinates via j(t, s(t, x), z(t, x)) = j(t, x) we obtain 

V a ,-j = div r „j + ^j-i/, (4.10) 

where divr ez j is the divergence of j on Y £Z . In the Appendix we compute 

A x b = A r J - i(« + sz\S\ 2 )d z b + ±d zz b + h.o.t. , (4.11) 

where k is the mean curvature (the sum of the principal curvatures) and |«S| is the 
spectral norm of the Weingarten map S. In addition we note that (see Appendix) 

Vr ez b(s,z) = V T b(s,z)+ h.o.t. , 
div r ez j(s,^) = div r j(s,z) + h.o.t. , 
A Te X s > z ) = A r b(s,z)+ h.o.t. 

where Vr,divr, Ar are the surface gradient, the surface divergence and the surface 
Laplacian on T. 
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4.3.2 Matching conditions 

We now assume an e-series approximation of the unknown functions (p, /i, v,p, w, . . . 
which in the inner variables we will denote by M, V, P, W, . . . . Denoting by $o + 
e$i + . . . the inner expansion and by ipo + sifi + ■ ■ ■ the outer expansion of the phase 
field we obtain the following matching conditions at x = 'j(s): 

lim = <A)(z±), (4.12) 

z— >±oo 

lim d z ^!(z,s) = V<f (x±)-u (4.13) 

z— >±oo 

where 09n(x±), . . . denotes the limit limo9n(^ ± 5u). In addition we obtain that if 

<5\0 

$1(2;, s) = A±(s) + B±(s)z + o(l) as z — > ±00 the identities 

A±(s) = ip x (x±), B ± {s) = Vy? (x±) • v (4.14) 

have to hold (see [17] . [2Tj). Of course similar relations hold for the other functions 
like v, )tt, 

4.3.3 The equations to leading order 

Plugging the asymptotic expansions into (I4.ip - (l4.5p we ask that each individual co- 
efficient of a power in e vanishes. The equation f)4.4p gives to leading order ~: 

= d zz % - V'($o) (4-15) 

and from f)4.12p we obtain 

$0(2) ±1 for z ^±00. (4.16) 
We now choose the unique solution of f)4.15p . (I4.16P which fulfills 

$o(0) = 0. 

We in particular obtain that $0 does not depend on t and s. Equation (14. 2p gives to 
leading order 

d s V o -v = d z {V o -v)=0. (4.17) 
The matching condition requires that (V ■ v)(z) is bounded. Hence 

(v -u )(x+) = hm(V ■ v )(z) = lim (V ■ u )(z) = (v • u )(x-) . 

z—>oo z— >— 00 

This implies 

[v -u]± = a, 

where [w]t(x) = u(x+) — u(x—) denotes the jump of a quantity at the interface. 
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The analysis of (14.31) now depends on the ansatz for the mobility. We have to 
distinguish between four cases. 

Case I : m £ (ip) = m . 

At order ^ we obtain from ( I4.3P 

= d z (m d z M u) ■ v = m d zz M . 

Matching implies that Mo is bounded and hence Mo is constant. In particular we 
derive 

[Hi = o. 

Case II : m e ({p) = em . 

At order | we conclude from ( I4.3P 

- V8 Z % + (v • v)d z % = d z (m d z M u) ■ u = m d zz M . (4.18) 

We obtain from matching 

d z M — > for z — > ±oo . 
Integrating (14. 18j) with respect to z gives 

V = v ■ v . 

In addition we get d zz M = and hence M does not depend on z. 

Case III : m e (<p) = - <^ 2 )+ . 
At order 4y we get 

= d z ( mi (l - <S> 2 Q )d z M u) ■ v = d z {m x {\ - &*)d z M ) . 

Matching implies 

mi(l - $l)d z M for z ->• ±oo . 

We hence obtain 

7^(1-^)9^0 = 0, 

which implies 

M = M (s,t). 

Case IV : m e {ip) = mi(l — y2 2 )+ . 
At order | we get 

-V<9 Z $ + (V ) • ud z $ = d z ( mi (l - ®l)d z M ) 
and hence combining arguments from the Cases II and III above we obtain 

V = v • v and M = M (s, t) . 
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We now analyze the diffusion equation for the soluble species. Equation (14. 5p gives 
to leading order \ 

d z (K(%)W d z (g\W + /3($ ))) = . (4.19) 

Matching to the outer solution gives that g'(Wo) + /3($o) is bounded. Multiplying 
(I4.19P by g'(Wo) + /3($o), integration and integration by parts gives 

K(^ )Wo\d z (g\W ) + (3(®o))\ 2 dz = . 

Assuming Wo > we obtain that 

g'(W ) + (3(%) 

does not depend on z. Hence 

[g'(w o )+P(cpo)}t = 0. 
For the choice g{w) = wlogw — w we obtain 

logw (a;+) - \ogw (x-) + /3 2 - fa = 
which yields the Henry jump condition 

— — — = exp(/3 2 - fa) ■ 
w {x+) 

Applying (14. 9 p for each component we obtain 
V x v = i^V<8»i/ + V r „V, 

D x v = ii(^ V ® l/ + iy ®^ V ) + l( V r £Z V +(Vr £Z V) T ). 

Defining £(A) = |(A + A T ) for a quadratic matrix A we compute 

V x -(ri(<p)D x v) = ^d z {r } {^)£{d z Y®u))u + \d z {r ] {^)£{Vv ez V))t' 

+lV r „ ■ (!j(*)£{d M V <8> u)) + Vr £2 • (v(®)£ (Vr ez V)) 
= ^d z ( V (<Z>)£(d z V ® v)v) + '-dM^iVr^v) 

+ l £ Vr £Z ■ toWSfrV <8> u)) + Vr £2 ■ )5(V r „ V)) , 

where we used c^i/ = 0. We conclude from (I4.17P 

(i/®^v )i/ = (a 4 Vo-i/)i/ = o. 

The fact that <3>o does not depend on t and s and (14. 9p imply 

Vyj <S> Vy> = js(d z $ ) 2 (v ® f ) + id z $i^$ (f ® f) + h.o.t. 
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and since {SJ T v)u = we get 

eV ■ (Vy? ® V<^) = j?d z (d z <S> ) 2 v + i(«9 2 $o) 2 (V r ■ u)u + ^(d^d^u + h.o.t.. 

Since the chemical potential to leading order does not depend on z we observe that 
the term div(v (g> m(</?)V/i) gives no contribution to the order \. Hence we obtain 
at the order from the momentum equation 

ad z {d z ^fu + d z ( V ($ )d z V ) = . (4.20) 

Multiplying (I4.20p with u, taking d z u = and <9 2 V ■ i> = into account gives 

ad z ((d z $ ) 2 ) = . 

Hence ()4.2Up implies 

dM®o)d z Vo) = 0. (4.21) 

Matching implies that V (z) is bounded. This implies that (14.211) interpreted as an 
ODE in z has only solutions Vo which are constant in z. This implies after matching 

[v o ]i = 0. (4.22) 

4.3.4 The generalized Gibbs— Thomson equation 

The equation for the chemical potential gives to the order e° 

M = a^'($o)$ i - &d zz $! + &d z $ Q K + P'($o)W ■ (4.23) 

In order to be able to obtain a solution $x from (I4.23P a solvability condition has to 
hold. This solvability condition will yield the generalized Gibbs-Thomson equation. 
We multiply (14.231) with <9 2 $o, integrate with respect to z and obtain (using the facts 
that M and Vo do not depend on z): 

/oo 
(V>"($ o )0*$o$i - d zz ^d z %)dz 
-oo 

/oo /*oo 
(d z %) 2 dz + / d z /3{%)W dz. 
-oo J — oo 

Defining Co := /^(c^o) 2 dz we obtain after integration by parts, using the fact that 
d z §o(z), d zz §o(z) decay exponentially for \z\ —> oo, 

/oo 
d z (i(/($o) - <9 22 $ )$idz + vc k 
-oo 

pOO POO 

+ / d t (P($o)W Q + g(Wo))dz- / (P($ )+g\Wo))d s W dz. 



oo 
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Since /3($o) + 9'(W ) does not depend on z and since ip'(&o) — d zz Q = 0, see ( I4.15p . 
we obtain 

2/i = (tk+ [g(w ) - g'(w )wo]t 

where a := co<x. This identity has been derived for a general function g. In the case 
g{w) = wlogw — w we obtain 

2/x = ok- [w ]t ■ 

4.3.5 Interfacial flux balance in the sharp interface limit 

We now expand the equations (14. 3 p and (I4.5P further in order to obtain contributions 
of the diffusive fluxes in the interface. The result will depend on the choice of the 
mobility. In the cases II and IV the interfacial diffusive fluxes for <p are scaled such 
that they do not contribute to a limiting sharp interface problem. We hence consider 
only the cases I and III. 

Case I : m e ((p) = m . 

At order | we deduce from (I4.3P 

(-V + V • u)d z % = dzimodzMt) . (4.24) 

Matching gives d z M\ — > V/io • v for z — > ±oo. Integrating (I4.24p gives 

2(-V + v -v) = m [V// • H± • ( 4 - 2 5) 
Case III : m e (<p) = ^(1 - ip 2 ) + . 

We have up to order e° in the interfacial region (setting m(ip) = mi(l — if 2 )) 

V • (m(tp)Vti) = ±d z (m(%)d z M ) + 

\d z {m{<$> Q )d z M x + m'($o)$i9*M ) + ±V r • (m($ )«9 2 M ^) 

+d z (m(%)d z M 2 + m'($ )$i3,Mi) + V r • (m^oRMif) + V r • (m($ )V r M ) . 

Using d z Mo = we obtain from (14.31) at order 4j 

U = d z {m{<& )d z M 1 ). 

Matching gives 

m(§ )d z M 1 — ^ for z — > ±oo . 
Hence m(<3> )<9 2 Mi = and 

Mi = Mi(s,t) . 

At order | we obtain from (I4.3p . using d z M\ = and d z M = 0, 

- V<9 Z $ + (V • v)d z $ = d z (m($ )d z M 2 ) + V r • (m($o)V r M ) . (4.26) 
Matching gives 

m($o)<9 2 M 2 -»• for z -> ±oo . 
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Integrating (I4.26P gives 

2(-V + v ■ i/) = / Vr-(m($ )V r Mo)tte. 



Since m($ ) does not depend on s, we obtain Vr • (m($ )VrM ) = m($ )A r M . 
Altogether we deduced 

2(-V + v • i/) = mA r /io (4.27) 

where m = m($o)<^- 

Finally we deduce the flux balance for w at the interface in the limit e — > 0. 
Equation (14.51) gives to order | 

- Vd z W + (V • f)<9 z Wo = d z (K(^ )W d z (g"(W )W 1 + /3'($ )$i)) . (4.28) 
Matching gives 

d z (g"(W )W 1 + /3'($ )$i) -)• V(^K) + /9(±1)) ■ f 
for 2 ±oo. Integrating (I4.28P gives 

(V - v • u)[w ]t = -[Kw Vg'(w )]t ■ v. 
In the case g(w) = wlogw — w this identity reduces to 

(y-v -u)[w ]i: = -[KVw ]t-u. 

4.3.6 The momentum balance in the sharp interface limit 

We now want to analyze the momentum equation to the next order. The term 
V • (r)(<p)Dv) gives to the order -, 

Matching requires lim d z V~i(z) = Vv (i±)y and hence 

z— >±oo 

d z \i (g)u + V r V -»■ V x v for z -»■ ±oo . (4.29) 

In the cases II and IV the term div(v £g> m s ((p)Vfi) gives no contribution to order - 
and hence we obtain for the momentum equation at order |: 

- 9 z (p($o)Vo)V + 9 z (p($o)(Vo®V ))^ 
29 2 (r / ($ )^(^V 1 ® - 20 a fa($o)£(VrV o )i/) 

- a(9 2 $ ) 2 ^ + ad s (d z $id g & )v + d z P v = . 

Integrating with respect to z gives after matching and using (14.291) 

-[p v ]iV + [p v ]lv ■ u - 2[T)E(V x v )]tv -a(J (3 Z %) 2 dz\ ku + \p ]±u = . 
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We obtain since v • u = V 

- 2[r)Dv ]tv + \po]±v = okv . (4.30) 

Case I (m e ((p) = m ): 

We note that Mo does not depend on z and that at x — 7(s) 

dgMtU + V r M -»■ V/i (x±) for 2 -»■ ±00. 
For the momentum balance we hence get to the order ^ 

^(p($o)V )V + ^(p($ )(V ® V ))u - <9 Z (V <g> i[p ]l(^Mn/ + V r M ))v 

- 2d z ( V ($o)£(d z V 1 <g> - 2^(77($o)^(V r V )i/) 

- a(d x $ ) 2 KV + <r<^(<^$A$o)f + ^Po^ = • 
Integrating this equation with respect to z we obtain since 

/OO 
d z (V ® |[p ]i(^Mn/ + V r M ))efe = v <g> §[p ]l [V/i ]l 
-00 

the following interface condition 

Miv (v • 1/ - V) - |[po]lv [Vp ]- • f - [2r/£)v ]^f + \po]tu = ctkv . 
Using f)4.25p we can simplify the above identity to 

- [277-Dv ]ii/ + \po]- v = OKV ■ ( 4 - 31 ) 

We hence observe that in case I the sharp interface problem leads to the classical 
conditions (14.221) . (I4.3ip for v and p on the interface. 

Case III (m £ (ip) = ^(1 - <p 2 ) +y ): 
Using d z M = d z M\ = we observe that div(v <g> m £ (<p) Vp) gives to the order - the 
contributions 

mi9,((l - ($o) 2 )+^M 2 (V ® 1/)) • v + div r (m 1 (l - $g)+(V ® V r M )). 

Taking the flux terms into account we obtain from the momentum balance at order 
1. 

- ^(p($o)Vo)V + ^(p($o)(Vo®V ))i/ 

+ §Mtmi^((l - ($o) 2 )+^M 2 (V ® 1/)) • v 
+ |[p ]t div r (mi(l - $g)+(Vo ® V r M ))) 

2^(r 7 ($ )£(^V 1 ® - 2^(^($o)5(V r V )i/) 

- a(9 z $o) 2 ^ + ^(^^$0)^ + 5 Z P ^ = . 
Integrating the above identity with respect to z gives after matching 

[pov ]l(v • v - V) - i[p ]i^div r (v <g> V r po) - 2[7pDv ]li>' + \po]tv = okv . 
Using ( I4.27P we can simplify the above identity to 

- |[po]lm((V r po) • V r )v - 2[r)Dv Q ]±v + \po]±u = gkv . (4.32) 
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Remark 4.1 The interfacial stress balances (I4.30P and ( I4.32p have a normal and a 
tangential part. For (I4.30p we obtain the classical conditions 

— 2u ■ [i]Dvo]~tv + [po]l — ok an d T ■ [t7-Dvo]Ii^ = 

where the last identity holds for all r G M. d with r • i> — and hence states the 
continuity of the tangential interfacial stress. 

The interfacial stress balance (I4.32p has the normal part 

-\[p Q \tmu ■ ((Vr/io) • V r )v - 2u ■ [vDv ]tu + [p ]i = ok 

and the tangential part 

— |[po]-^T • ((Vr/Wo) • Vr)v — 2r • [i]Dv }^i/ = for all tangential r, 

i.e., in this case the stress [rjDv^v can be discontinuous also in tangential directions. 



5 Free energy inequalities for the sharp interface 
limit 

The sharp interface limit is now given as follows. We search for a smoothly evolving 
hypersurface V which for all t > separates Q into open sets f2_(t) and and 
we look for functions v,p,w which are all defined for x G £1 and t > 0. In the cases 
I and III we seek in addition a chemical potential /i which in case I is defined on 
Q x (0, oo) and for case III the potential /i is defined on the interface T only. 

Assuming for g the form g(w) = w(logw — 1) we need to solve in the cases II-IV 
the system 

Pi{dtv + div(v <g> v)) — r/jAv + Vp = 0, 

divv = 0, 
d t w + v • Vu> = KiAw 

in the bulk regions f2_ and In case I the momentum balance has to be replaced 
by 

Mv + div (v ® ( Pi v + e^moVfi)) - r^Av + Vp = (5.1) 
and in addition we need to solve 

A/i = 

in the bulk. On the interface T we require 

0, 

exp(/3i -fa), 
-[KVw]t ■ v 



(V 



[vjl 

Wi 
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where W\ = w(x—) and w 2 = w(x+). All other conditions depend on which of the 
cases I-IV we consider. 

Case I (m e (ip) = m ): In this case we require on the interface 

- [2T)Dv\tv + \p\±v = <tku, (5.2) 
2(vi/-V) = m [V/i]l-i/, (5.3) 
2 fx = ok- Hi . (5.4) 

Here, it is remarkable that the condition (15. 2p is not affected by the diffusional flux 
and that the interface is not passively transported by the interface, i.e. v • v ^ V is 
possible, compare Abels and Roger [3]. 

Cases II and IV : 

— [2i]Dv]t_ ■ v + [p}^_u = okv , 
V = Vi/, 

i.e., in this case we recover a standard free boundary problem for the Navier-Stokes 
system which in addition is coupled to the flow of a soluble species. 

Caselll (m e { V ) = ^{l- V 2 ) + ): 
We require that 

-||>]im((Vr^)-V r )v-2[77Dv]ii/+[p]ii/ = anv , (5.5) 

2(vi/-V) = mAr/i, (5.6) 

2fi = ok - [w]± . (5.7) 

In this case the diffusion of the two components is limited to the interfacial region. 
In fact the well-known surface diffusion flow V = — A r K, see [29] |P3 EH ES], is in the 
new model coupled to fluid flow. In this case we observe that the surface flux enters 
the interfacial stress balance. 

In all cases we require v = 0, Vw • u — on dfl and in case I we in addition 
require V/i ■ v = on dQ. In order to derive a free energy inequality for the sharp 
interface limit we need the following transport identities, for a proof see [TU] . 

Lemma 5.1 (Transport identities). Let T = (T(t)) t >o with T(t) C £1, for all t > 0, 
be a smooth evolving hypersurface and let f be a quantity which is smooth for all 
t > and x G Q \ T(t) and such that [/]! exists at T. Then it holds 

^- I lds x = - [ kV ds x = - [ {kv) ■ (Vu) ds x 
dt Jr(t) Jr(t) Jv(t) 

and 

dt 



— / fdx 



[ d t fdx- [ [f} + Vds x 
Jn Jr(t) 
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We are now in a position to compute the dissipation rate and in conclusion derive 
a free energy inequality. 

Theorem 5.2 (Free energy inequality) . 

A sufficiently smooth solution of the sharp interface problem with T(t) C Q, for 
all t > 0, fulfills 



d 
~dt 



I (ttI v | 2 + (9(w) + f3w)) dx+ / ads 
Jn^^ ' Jr(t) 



= -V < , 



provided the integrals are finite. Here p = p±, p 2 and (3 = (3i, (3 2 in the two phases. 
The quantity V is given as 

Case I : V = I 2r]\Dv\ 2 dx + / K±\ Vw\ 2 dx + / m \Vp\ 2 dx, 

Jn Jn w Jn 

Case II and IV : V= [ 2r]\Dv\ 2 dx + I K±\Vw\ 2 dx, 

Jn Jn w 

: V= 2r]\Dv\ 2 dx+ / K±\Vw\ 2 dx+ I m\V r p\ 2 ds x . 
Jn Jn w Jr(t) 



Case III 



Proof: Let us first discuss the terms which have to be treated in the same way 
in all four cases. The total free energy of the soluble species fulfills (using g{w) = 
w log w — w) 

"7" / (g(w) + /3w)dx = [ d t w(\ogw + (3) dx — f [g(w) + f3w]tVds x . 
dt Jn Jn Jv(t) 

Since d t w + (v • V)u> = KV ■ (toV logw), we obtain 

d t w(logw + (3)dx = / (logw + /3)(-(v • V)w + KV ■ (wV\ogw))dx 
Jn Jn 

— — I (v • V)(g(w) + f3w) dx — I Kw\ V log w\ 2 dx 

Jn Jn 

— / [(w log w + (5w)KV \ogw\Z ■ v ds x 
Jv{t) 

— — I Kw\V \ogw\ 2 dx + I [g{w) + (3w]Zv ■ v ds x 

Jn Jr<t) 



[(\ogw + (3)KVw]t -vds x 



r(t) 



— — KwlV logw\ 2 dx + / [g(w) + f3w]Zv ■ vds^ 
Jn Jr(t) 

+ / (logw + p)(V-v-v)[w]zds x 
Jth) 
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where we used divv = and the fact that Henry's law implies the continuity of 
logw + p. In addition we have 

— f ads x = — [ anVds x . 



dt J T 



We now give a proof for the cases I and III in detail and discuss the cases II and 
IV afterwards. 

Case I: We introduce the abbreviation J = — p2 ~ pl m(<p)Vp. For the kinetic 
energy we compute, using a transport identity and ( 15. ip 



d 
dt 



I p^-dx= ( p\-d t vdx- I [p}t^-Vds x 
Jn 1 Jn Jr(t) * 

= / v • (— div(v <g> (pv + J) + 2p div Dv — Vp) dx — f [p)t^-V ds a 
Jn Jr<t) 2 



Using div v = and div J = in £l± we obtain 

v ■ div(v ® (pv + J)) = |V|v| 2 ■ (pv + J) 
Integration by parts on Q + and now gives 



d f |v 

dl 



f p^Ldx = - [ 2 V \Dv\ 2 dx + [ I| v | 2 ([p] + (v v- V) + [J]l ■ u)ds x 
Jn 2 J n J r ^ 

+ v • (\p]±u - [2r]Dv]tiy)ds x 

Mt) 



= — 2rj\Dv\ dx + aK~v-vds x . 
Jn Jr(t) 

where we used the stress balance (I5.ip . the interface conditions (15.21) . (15. 3p and the 
fact that v is continuous. 

We compute using integration by parts, (15.31) and Ap = in the bulk 



/ 2(v ■ u — V)pds x = / m p[Vp]^ • v ds x (5.8) 

M*) Mn 

= — V • (pm Vp) dx = — l m \Vp\ 2 dx. 
Jn Jn 



Altogether we obtain using the generalized Gibbs-Thomson law (15.41) and the 
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fact that g(w) + (3w — (log w + (3)w = —w 

J (p^~Y~ + d( w ) + Pu^j dx + J vds^j 

- I 2rj\Dv\ 2 dx - / K±\Vw\ 2 dx 
Jn Jn w 



d 
It 



(5.9) 



+ / (v-u-V)(aK-[w)z)ds a 
'r(t) 



= - f 2r)\Dv\ 2 dx- I K±\Vw\ 2 dx- [ m \Vfi\ 2 dx. 
Jn Jn w Jn 

Case III: Using ( 15. 5p . ( I5.6p . ( 15. 7p . setting Jr = pl ~ p2 mVryU and using the divergence 
theorem on manifolds we compute 



- [ 2r]\Dv\ 2 dx+ [ ||v| 2 [p]l(v v- V)ds 2 
Jn Jr(t) 

+ / v ■ (\p]+u - [2rjDv]±v)ds x 
Jrtt) 



2r]\Dv\ 2 dx 



+ / v ■ (\p)±v - [2r]Dw]tf - |vdiv r J r )ds a 

r(t) 



2rj\Dv\Mx+ I v {\p\Zv- [2r/Dv]+^ + (J r ■ V r )v)ds 3 
r(*) 



/ 2ri\D\\ 2 dx + / <7K\-vds x . 
Jn Jr{t) 



We now obtain 

d 
dt 



(^-\v\ 2 + g(w) + (3w) dx + a ds a 
V2 / J v{ t) 

- / 2ri\Dv\ 2 dx - / K±\Vw\ 2 dx 
Jn Jn 



+ (vi/-V)(<TK-[w)z)ds 3 
'r(t) 



— — I 2r]\D\\ 2 dx — / iC— | Vu>| 2 + / m(Arn)/J>ds x . 
Jn Jn w Jr(t) 

Applying the divergence theorem on manifolds for the last term now gives the result 
in case III. 

In the cases II and IV we have V = v • i> and the calculations simplify. In par- 
ticular we do not need to treat the term J^/ t \(v ■ v — V)2fids x and diffusion plays no 
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role for the fluid flow. 



Remark 5.3 Using formal asymptotic expansions one can also derive a sharp inter- 
face limit for the non frame invariant model f l2.38p - fl2.4ip . In this case the flux term 
involving the gradient of the chemical potential V/x does not enter the momentum 
balance law ( 15. ip . The interface law (I5.2p has to be replaced by 

[p]^v(v • u — V) — [27]Dv\^h> + \p]~tv = OKV 

and the summand involving Vp£t in (15. 5p drops out. In addition the generalized 
Gibbs-Thomson laws f !5.4p and (15. 7p now contain an additional term involving the 
velocity and are given as follows 

2fi = ctk- ||v| 2 [p]+ - [w]t- 

We also remark that the free energy inequality in Theorem 5.2 also holds in this case. 
We refer to [3] for details. 



Appendix 

We use the notation of Section H~31 and prove the identities (14.91) . 14. lUp and f)4.1ip . 
Let (si, . . . , Sd-i) £ U and G £ x (s, z) = 7(5) + ezu(s), where here and in what follows 
we omit the t-dependence. Then 

<9 S1 7 + ezd si u, d Sd _ x ~f + ezd Sd _ x v, ev 

is a basis of M d locally around T. We define the metric tensor in the new coordinates 
as follows 



9ij = (d Si j + ezd Si u) ■ (d Sj j + ezd Sj u) i,j = l,...,d-l 
gid = g<u = (darf + ezd Si v) ev = % = 1, . . . , d - 1 , 



9dd 



ev ■ eu = e 



where we used that d Si v ■ v = ^d Si \is\ 2 = 0. We set 



Gez = (gij)i,j=l,...,d, Gez = (9ij)i,j=l,...,d-l, (Gez) = (g V )i,j=l,...,d, (Gez) = (g l3 )i,j=i 



4-1 



and hence 



\ 



\ • ■ • e 2 J 



( 



\ 



, (Gs 



,-1 



{GJ' X 



\ 0---0 e~ 2 J 



35 



Denoting by Sd the ^-variable we have for a scalar function b(t, x) = b(t, s(t, x), z(t, x)) 

d d-l 
i,j=l i,j=l 

= VrJ + l £ d z bu 

where we used the fact that = for i = 1, . . . , d — 1. Here Vr £Z & is the surface 
gradient Vr ez &|r ez on T ez := {"f(s) + ezu(s) \ s 6 U}. In addition we compute for a 
vector quantity x) = s(t, x), z(t, x)) 

d d-l 

div.j = ^^aj-5. J C^=^^'aj.a. j G; + ^j.ei/ (5.10) 

ij'=l ij'=l 

= div r £Z j + \d z ] ■ v , 

where divr ez j is the divergence on T £Z . We remark that Vr ez &- v = , as u is normal 
to Y £Z . We hence obtain 

d z {V T J-v) = 

and 

^(V r „6)-i/ + Vr„6-9*i/ = 0. 

Since <9 2 f = 0, we get 

(<9 z V r J)-i/ = 0. 

We now compute 

A,6 = div x (V x b) = div x (V T J + \d z b v) 

= div ra ,(Vr„S) + KVr.A6)'«' 

+j9*6V r „ • ^ + J) • " + • + jsd z bd z v ■ i/ . 

Because of (Vr ez <9 2 6) • i/ = 0, (<9 z Vr £Z &) ■ u — 0, u — V x d, ( I5.10P and <9 z i/ = 0, we 
obtain 

A^fe = A r J + i^d^S + ±d zz b . 
Introducing g £z = det Q £Z we get 

d-l 

Since 

<fa = 9 Si 7 • 9 Si 7 + h.o.t. , i,j = 1, . . . , d - 1, 
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we derive 



Vr E j(s,z) = V r b(s,z) + h.o.t. , 
divr ez j(s,z) = div r j(s,z) + h.o.t. , 
A r£z b(s,z) = A r S(s, z) + h.o.t. 

where Vr, divr and Ap are computed on T EZ with the metric tensor Q . From Gilbarg 
and Trudinger [22J (Lemma 14.17) we get (denoting by Ki the principal curvatures) 

d-l 

'^^6K^Z + h.O.t. 

i=l 

where k is the mean curvature and \S\ is the spectral norm of the Weingarten map 5. 
Hence we obtain 

A x b = A r b - i(« + £z|S| 2 )<9 2 6 + ^«9 22 6 + h.o.t. . 



d-l d-l d-l 

i=l i=l i=l 

= —K — ez\S\ 2 + h.o.t., 
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